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ABSTRACT 

Central difference approximations to the fluid dynamic equations require an artificial 
viscosity in order to converge to a steady state. This artificial viscosity serves two purposes. 
One is to suppress high frequency noise which is not damped by the central differences. The 
second purpose is to introduce an entropy-like condition so that shocks can be captured. 
These viscosities need a coefficient to measure the amount of viscosity to be added. In the 
standard scheme, a scalar coefficient is used based oa the spectral radius of the Jacobian of 
the convective flux. However, this can add too much viscosity to the slower waves. Hence, 
we suggest using a matrix viscosity. This gives an appropriate viscosity for each wave 
component. With this matrix valued coefficient, the central difference scheme becomes 
closer to upwind biased methods. 


*This research was supported by the National Aeronautics and Space Administration under NASA Con- 
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Science and Engineering (ICASE), NASA Langley Research (tenter, Hampton, VA 23665. 
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I. Introduction 


In recent years, central difference schemes have been used with much success to solve 
transonic flow problems about aerodynamic shapes. These schemes are second order ac- 
curate for sufficiently smooth meshes and have an added artificial viscosity to stabilize 
the scheme and reach a steady state. This artificial viscosity is usually a blend of two 
terms. One is a fourth difference that stabilizes the even-odd modes that appear with 
central differences and constants coefficients. Without this viscosity, one cannot reduce 
the residual beyond some level because of a remaining high frequency mode. The second 
viscosity term is a nonlinear second difference that limits oscillations in the neighborhood 
of shocks. A nonlinear shock detector preserves the second order accuracy of the scheme 
in smooth regions. 

An advantage of the artificial viscosity approach is that it allows the user control 
over the amount of dissipation. Nevertheless, one sometimes finds that there is too much 
dissipation in the numerical solution. Changing the global constants that appear in the 
formulas is not sufficient to construct an artificial viscosity that is appropriate in both the 
shocked and smooth regions of the flow. For some problems, we need to severely limit 
the viscosity in some smooth regions, e.g., near the trailing edge, while still maintaining 
stability near the shocks. Hence, although the standard artificial viscosity works well in 
most cases, it is not sufficiently flexible to handle more delicate problems. 

In order to improve the existing artificial viscosity, we shall make use of ideas used in 
the construction of upwind schemes. In particular, we shall replace the scalar coefficient in 
the artificial viscosity by a matrix. To prevent difficulties near stagnation points, a cutoff 
is introduced that depends on the spectral radius of the matrix. By varying the cutoff, 
one can obtain an appropriate average between the original scalar viscosity and the new 
matrix viscosity. Since the matrix viscosity reduces the amount of smoothing on the slower 
waves, it will improve the total accuracy of the scheme. 


II. Finite Volume Formulation and Artificial Viscosity 


The Euler equations for an inviscid compresiible flow can be written in divergence form 


as 


dt dx dy dz 


( 1 ) 


where 


Q = {p,pu,pv,pw,EY 


(2a) 



/ = (p«, P« 2 + p , puv, puw, ( E + p)u)‘ (26) 

g = (pv, puv, pv 2 + p, pvw, ( E + p)v)‘ (2c) 

* = {pw, puw, pvw, pw 2 + p.^ + pjtv)' [2d) 

and for an ideal gas 

P= (7-l)[^-p(u 2 + v 2 + tv 2 )/2]. (2e) 

We can also write (1) in the form 

dQ 

~m + d,v(F) = °- (ifr) 

We integrate (l) over a three dimensional cell and consider Qx,j,k <is an approximation 
to the average of Q over the cell. Hence, 


dQij.k fff divFdV 
dt fffdV 


or using the divergence theorem, 


+ J j & • ndS = 0 . ( 3 ) 

Hence, the time change of the average Q is governed by the fluxes entering and leaving the 
cell. 

This finite volume approach leads to a pure central difference method for Cartesian 
grids. Though this scheme is stable for constant coefficient hyperbolic equations it is 
subject to instabilities that will prevent the convergence to a steady state. To enhance 
this convergence a fourth difference viscosity is added to the scheme. The fourth difference 
causes oscillations in the neighborhood of shocks. Hence, a shock detector is constructed 
and near the shocks the fourth difference is turned off and only the nonlinear second 
difference is operative. The total artificial viscosity, V , is the sum of such second and 
fourth differences in each coordinate direction. 


tot 









( 4 ) 


Hence it is sufficient to describe these terms in the £ direction. Since we only take differ- 
ences at neighboring points the artificial viscosity is always in conservation form. 

The first difference is defined as 


Qi+l ,},k Qi,i,k 


(5a) 
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and the second £ difference is defined as 


Ei.j.k — Di+lj'k Di_i j k . 


(56) 


We then form the second and fourth differences. In particular the fourth difference is 
formed as a second difference of a second difference with positive weights [3,8]. Hence, 




D 






( 6 ) 


Let, 


Pi+l,j,k ~ 2 Pi j,k + Pi-1, j,k 
Pi+l,j,k + 2Pi j,k + Pi-l,j,k 


(7a) 


v: i k is used to detect the location of shocks. When v iijik is large then the fourth difference 


is reduced. Let, 




= If (2) max(K-i i)it) v i jt k, Vi+i'j'k, Vi+ 2 ,j,k)- 


( 8 ) 


We also multiply a by a function of the Mach number to reduce a near the surface. Let, 
A = = |^, where F, G,Tt are the fluxes in the coordinate system (£,r?,f). 

Let A be a measure of the fluxes. The original code chose A as 


A« = A" = A f = p(A) f p(B) + p(C) 


(9a) 


where p is the the spectral radius of the matrix. For problems with a highly stretched 
mesh it was found [2,3,8] that for increased accuracy one should choose 

A< = p(A), A" = p(l 0, A ? = p(C). (96) 


are constants that determine the level of the second and fourth differences. 
These constants are given as input to the code. Then 



£, + "^ , + J . ',k^i+ \,j,k 

= A tiii *max(0, K (4) - o i>j<k ). 


(10a) 

(106) 


In order to imitate the upwind type [7] algorithms we now replace the scalars in (9b) 


by matrices. Hence, 

A« = |A|, A" = |JJ|, A f = |C| (9c) 


where \A\ = r|Af|T _1 when A = TA^T~ l and Ac is a diagonal matrix with the eigenvalues 
of A as its entries. This definition of A can lead to difficulties when an eigenvalue is near 
zero. Hence, we modified A € to be 


A^ = dta< 7 (max(aj,< 7 />(.A)) a, = e.v. of A 


( 11 ) 


3 


where q is a specified constant. When q = 1, then A f = p(A) • I and so (9b) is recovered. 

When 9 = 0 then no modification of A^ is done. In general, we found that q = 0.2 gives 
good results. 

We point out that the use of (9c) does not allow for a constant enthalpy solution and 
so enthalpy damping cannot be used [5]. 

III. Results 

We consider the central difference code with Runge-Kutta time stepping [3,5]. As de- 
scribed above we use a matrix valued artificial viscosity which approximates TVD type 
schemes [6,7,9]. The fourth difference viscosity is still needed to allow the multigrid accel- 
eration to quickly reach a steady state. We consider inviscid flow about a NACA0012 with 
M °° = 0.8, a = 1.25°. A 192 x 32 C mesh is used with 128 points on the airfoil. In [l] it is 
pointed out that the standard code smears the weak shock on the lower surface. In Figure 

1, we plot the result for the standard scheme, but without enthalpy damping. In Figure 

2, we show the same case but using the matrix viscosity. The convergence rate is slowed 
down since the fourth order viscosity is not as strong but the shock on the lower surface is 
sharper. There is an overshoot on the shock on the upper surface. This is due to fact that 
the cutoff (10b) is not sufficiently sharp. One way to improve this is to replace (7) by 

U i jk = |P»+l.J,t ~ Pi,j,k\ ~ \Pi,j,k ~ Pi-i,j,k\ . . 

|P»+l,/,fc ~ Pi,j,k\ + \Pi,j,k ~ Pi—l,j,k | + £ 

so that Uij k is one at discontinuities. One can also use the matrix coefficient only for 
the second difference but use a scalar coefficient for the fourth difference viscosity. This 
accelerates the convergence to a steady state but smears the shocks. The results presented 
used a four step Runge-Kutta algorithm with the artificial viscosity frozen after the first 

stage using a matrix viscosity rather than a scalar viscosity adds about 60% to the total 
CPU time. 
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Figure 2: Matrix viscosity 



Appendix 


We present the matrix viscosity in explicit form for three dimensions in general curvi- 
linear form. Let 

t=t{z,y,z), rj = rj{x, y, z), <T = f(x,y,z) 

be a change of variables from physical space (x, y, z) to a computational space (^, t?, c) such 
that the curvilinear mesh is mapped to a cube. Then (lb) can be rewritten as 

dQ dF dG dH 

dt + ( A1 ) 

where 

F = U + t,g + t.h 

G = V xf + r] y g + rf z h (A2) 

^ $zf d* $u9 ~t~ 

We next express (Al) in quasilinear form 


Let 




A 


— b= — r = ^L 

dQ dQ 8Q ’ 


(A3) 


9 = + i y v + £ z u> 

T = 9 x« + r] v v + T) Z W 

s — + CyV + $zU> 

be the three contravariant velocities. Also define h = ^ as the total enthalpy and let 
V 2 = te±l(u 2 + v 2 + w 2 ) then 


( 0 


ai 


o 2 


a 3 


0 ^ 


A = 


diV 2 — uq 
0.2V 2 — vq 
a 3 V 2 — t uq 
k -q{h-V 2 ) 


9 ~ (7 - 2)aju 
diV — (7 — l)a 2 ti 

diw — (7 — l)o 3 u 
d\ h — (7 — 1 )qh 


a 2 u — (7 — l)a!W 
9 ~ (7 “ 2)o 2 v 
a 2 w - (7 - l)a 3 t; 
d 2 h — (7 — l)gv 


a 3 u - (7 - l)aity (7 — 1)0! 
a 3 v - (7 - l)? 2 to (7 — l)a 2 
9 -(7 -2 )a 3 w (7-1)03 
a 3 h — (7 — l)qw 
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19 J 
(A 4) 


For B we get a similar matrix with = »?*, ^2 = »7y> a 3 = 'l* an< * 9 replaced by r while 

for C we have a, = a 2 = fr, a S = £* and q replaced by s. 

Hence, we can find the absolute value of all three matrices in the same way. Let us 

assume that A has eigenvalues 

Ai 

A 2 0 

As 

0 As 

As 




where A x = q + yja\ + a\ + afc, A 2 = q - y[a\ + a\ + a\c and A 3 - q. In practice Ay, 

redefined by (11) to prevent A i from approaching zero. 

In order to find \A\ it is easier to use a two step procedure. Let 


are 


A\ = TyATi 1 


(A6) 


be a symmetric matrix. Since we can symmetrize A,B,C simultaneously the same T x will 
work for all three matrices. 

( Yl —(-7 — l)u/c -('y-l )vfc -{i-tfw/c (7 - l)/c \ 


Ti = -v 


-c 2 + 



0 

0 

0 

1 

w 



9 


As expected Ti, 7\ 1 do not depend on the matrice 

Oi, 

tt 2 , CL — 


( 

q ai c 

a 2 c 

a 3 c 

o \ 



9 

0 

0 

0 

A x = 

a 2 c 

0 

9 

0 

0 


a 3 c 

0 

0 

9 

0 



0 

0 

0 

0 

9 / 

Since A x (and B u C x ) is symmetric we can diagonalize it witl 

depend explicitly on a x , a 2 , a 3 and so is different for A u . 

B 1 ,C 1 

then 









A 

a 1 

g 2 

ds 

0 



-A 

a 1 

fl -2 

a 3 

0 

1 







■s 

II 

r« 

0 

Xi 

x 2 

X 3 

0 



0 

Vi 

y 2 

ys 

0 



^ 0 

0 

0 

0 

y/2A 


(AS) 


(A 9) 


and 


T; 1 = T. I 


The (xj,yj) are numbers that satisfy the following equations 

x\ + x\ + xl = y\ + y| + y\ = 2 A 2 
d\Xi 4 - 0,2X2 + 03X3 = UiJ/i + U2J/2 "h U3J/3 = 0 

^ll/l + Z2J/2 + X3J/3 = 0 

Zi!/i + 2^2 + x 3 y 3 = 0 
X\X 3 + I/lJ/2 = —20i02 
X\X 3 + t/ij/3 = —2a x a 3 


(A10) 


^2^3 + 2/2 J/3 — — 2a 2 as. 

It is difficult to give explicit formula for the in all cases since some of the a ; may be 

zero as long as a\ + a\-\- a\ 0. However, the final formula does not depend on explicitly 
knowing the Given T2 we find that A (A5) is given by 

A = T 2 A 1 T 2 ~ 1 = T 2 T 1 A{T 2 T 1 )- 1 . (All) 
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We now reverse the process and define 


( l*i 


w = 


|a 2 | o 
I'M 

0 | A3I 


IM ) 


where A ,• can be modified eigenvalues of A. Then 

\A\ = (T 2 T l )~‘MT 2 T 1 ). 


Let 


<Ti = 


Ai A 2 Ai A 2 


< r 2 = 


and define the row vectors 

«i = (7-l)( 


U 2 + V 2 + w 2 


,-u,-v,-u>, 1) 


We then have the matrices 


£2 = (-g, a 1,03,03,0). 


7 = 5x5 identity member 

Zi = {l u ul u vt u wl u hiiY 
Zi — (£ 2 ,u£ 2 + alii, v £ 2 + a 2 £i, + fl3£i, h£j + ?£i) t 

Z 3 — (0, ai£ 2 ,o 2 £ 2 ,a3£ 2 ,9£ 2 ) < 

and finally, 

w = As/ + (^V^) Zi + + Zs 

Because of the simple nature of the matrices Zj it is easy to multiply | A| times 
Define, 

ri = (l,tt,u,«;,A) < 
r 2 = (0,ai,a 2 ,a 3 ,g) 4 . 

Let ( , ) denote the standard inner product, then 

(£i,rx)=c 2 

(£i> r 2) = 0 


(A12) 

(A 13) 
(A14) 

(>115) 


(A17) 
a vector. 
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and if x is any column vector, then 


(^2> r i) = 0 
(t 2 ,r 2 ) = yp 

\A\x = X 3 x+ [^ " 1 c2 (*!,*) + r l+ + - gl ^ 2 ~ ( £ 2> X ) f 2- U 18 ) 
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